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The propagation of ion acoustic shocks in nonthermal plasmas is investigated, both analytically 
and numerically. An unmagnetized collisionless electron-ion plasma is considered, featuring a su- 
perthermal (non-Maxwellian) electron distribution, which is modeled by a n- (kappa) distribution 
function. Adopting a multiscale approach, it is shown that the dynamics of low-amplitude shocks 
is modeled by a hybrid Korteweg-de Vries - Burgers (KdVB) equation, in which the nonlinear and 
dispersion coefficients are functions of the n parameter, while the dissipative coefficient is a linear 
function of the ion viscosity. All relevant shock parameters are shown to depend on n: higher de- 
viations from a pure Maxwellian behavior induce shocks which are narrower, faster and of larger 
amplitude. The stability profile of the kink-shaped solutions of the KdVB equation against external 
perturbations is investigated and the spatial profile of the shocks is found to depend upon and the 
role of the interplay between dispersion and dissipation is elucidated. 



I. INTRODUCTION 

The study of shock waves in collisionless plasmas 
[l[ has received, both from the experimental 0-0| and 
theoretical |8l-ll4| point of view, a great deal of atten- 
tion in the last few decades. This widespread interest 
is justified by the central role that such structures 
play in many significant astrophysical scenarios such 
as cosmic rays generation [IHHll during supernova ex- 
plosionsJl7| , strong fields generation in the bow-shock 
region [18| and nonlinear dynamics of the solar wind 
[ljj, just to cite a few. Collisionless shock waves are 
routinely detected also in laser-plasma experiments 
which have therefore been proposed as smaller scale 
test-beds to systematically study astrophysical scenar- 
ios in the laboratory [2(| ■ The increasing level of detail 
offered by both laboratory measurements and space 
observations suggest the need for continuous refining 
of the underlying theory. Most of the theoretical stud- 
ies reported in the literature exploit the Korteweg-de 
Vries/Burgers (KdVB hereafter) equation which ele- 
gantly combines an interplay among nonlinearity, dis- 
persion and dissipation effects in the generation and 
evolution of shock structures. However, the outstand- 
ing variety of plasma scenarios in which shocks may 
occur has made it difficult to develop a unified theo- 
retical description of these intrinsically nonlinear ex- 
citations. 

A major theoretical challenge in plasma model- 
ing arises from the fact that most (e.g., astrophysi- 
cal) plasmas may present significant deviation from 
a Maxwellian behavior, mostly due to the presence 
of accelerated electron populations [2l| - |25| . Non- 
Maxwellian plasma particle distribution also occurs 
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in the lab [26l - [28| . e.g., during high- intensity laser- 
matter interactions, where the generated plasma 
clearly fails to thermalize within the short observa- 
tion time interval, a fact reportedly affecting the mea- 
surable characteristics of nonlinear structures propa- 
gating in the plasma [13, HH . Fundamentally speak- 
ing, such energetic particles populate the superther- 
mal region of the velocity distribution, and thus result 
in a long-tailed distribution function. A power-law 
behavior thus arises, at high velocities, a fact which 
motivated Vasyliunas [29| to introduce a phenomeno- 
logical long-tailed distribution, termed the kappa (k) 
distribution (after the real valued n spectral parame- 
ter involved in it), which succeeds in reproducing the 
observed power-law dependence at high energies [2!| , 
while remaining close to the Maxwellian curve at low 
(subthermal) velocities. In the years that followed, 
the kappa distribution has been proven remarkably 
successful in fitting space observations (see, e.g., the 
various references in Refs. |30M32|) and, more recently, 
laboratory measurements [26l428l |. Understandably, 
kappa theory has inspired a number of theoretical 
studies [3314371 ] , with the ambition of elucidating the 
effect of superthermal populations on the plasma dy- 
namics. 

Despite its success on the space-observational test- 
bench, Vasyliunas' apparently ubiquitous distribution 
function remained in the grey area of phenomenology, 
lacking a rigorous foundation from first principle. In 
the meantime, alternative forms of the original kappa 
function have appeared [HI, HH, citing all of which 
clearly goes beyond our scope here. The interested 
reader is referred to the interesting review and ex- 
tensive discussion presented in Ref. [3l|. A challeng- 
ing non-Maxwcllian-Boltzmannean paradigm appears 
to come from a different perspective, in Tsallis' non- 
extensive statistical- mechanical framework [3!l H(| , a 
recently proposed generalization of Boltzmann statis- 
tics, which seems to incorporate a parametrized fam- 
ily of non-Maxwellian functions which correspond sta- 
tionary equilibria which in fact maximize a generalized 
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entropy, in analogy with the Boltzmann function and 
Gibbs entropy (smoothly recovered from the Tsallis 
theory at some limit). Not surprisingly, the Tsallis 
theoretical club has started to attract its own mem- 
bers within the plasma physics community [U lU |42j • 
As a natural consequence, the relation between kappa 
and Tsallis theories was investigated in Space plasmas 
[43l [H| and, in fact, it was recently argued [3l| on a 
rigorous basis that the very foundation of the kappa 
distribution lies in the Tsallis theory. Although this 
remains an open topic in the community, it would ap- 
pear that the ubiquity of kappa distributions implies 
its deep foundations in Nature, which are still to be 
traced. 

Our target in this article is to investigate, from first 
principle, the effect of excess superthermality on a 
fundamental problem in plasma physics, namely the 
propagation of electrostatic shocks evolving on the 
ionic scale. Ado pting the re-theoretical framework in 
its original form [23 . |34| and relying on the analytical 
toolbox presented in Ref. H3 (details of which are to be 
omitted) , we shall develop a comprehensive formalism 
for one-dimensional ion-acoustic shock waves prop- 
agating in unmagnetized electron-ion plasma. Dis- 
sipation effects will be included by assuming non- 
negligible ion viscosity. The electron nonthcrmality 
is modeled by adopting a re-distribution for the elec- 
trons. A nonlinear evolution equation for the elec- 
tric potential will be derived, by using a multiscale 
(variable stretching) iterative scheme. The competing 
nonlinear and dispersive term(s) will be shown to be 
significantly dependent upon the re parameter which 
therefore significantly modifies the shock amplitude 
and width (compared to the Maxwellian case). An- 
alytical and numerical techniques will be employed 
to solve the associated KdVB equation (discussion 
above), and particular attention will be paid to the 
stability profile and the geometric characteristics of 
the shock excitations. 

The layout of this article is as follows. The build- 
ing blocks of the analytical model are presented in 
the following Section [XT] and its dispersion properties 
are analyzed, in Section IIIII A multiscale perturba- 
tion technique is then employed in Section IIV1 lead- 
ing to the working horse of the study, a Korteweg- 
de Vries/Burgers partial-differential equation (PDE) 
for the electrostatic potential. The shock-profile solu- 
tions of the K-dVB equation are presented and their 
dynamics is investigated in Section [V] Our results are 
finally summarized in the concluding Section IVII 



II. THE MODEL 

As mentioned in the introduction, we consider an 
unmagnetized electron-ion plasma where the electrons 
are assumed to be described by a re— distribution func- 
tion. The ions are assumed as a cold (Tj = 0) viscous 
population whereas the electrons will be considered 



as a fast neutralizing background. The electron iner- 
tia can in fact be neglected if we consider that ion- 
acoustic electrostatic waves move at phase velocity 
(v P h) which is much higher than the ion thermal speed 
and yet in turn much lower than the electron thermal 
speed: Vth,i *C v p h <C Vth,e- In a one-dimensional fluid 
description, the dynamics of the ions is governed by 
the following (normalized, dimensionless) set of fluid 
equations: 
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duj 
dt 
d 2 4> 

dx 2 



dx 
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(1) 
(2) 
(3) 



Here rij, Uj, (f> and 77 represent the ion density (nor- 
malized to the equilibrium ion density), the ion fluid 
velocity (normalized to the ion sound speed c s = 
(tcsTe/mi) 1 / 2 ), the electrostatic potential (normal- 
ized by ksTe/e) and the ion kinematic viscosity, re- 
spectively. Time and length are normalized in terms 
of the inverse of the ion plasma frequency (uj^ 1 = 
[47re 2 nio/ni.;] -1 / 2 ) and the plasma screening length 
(Xd = [ksT e / A-Ke 2 n e0 ] 1 / 2 ) respectively. We note 
that n e Q = riio was assumed to ascertain neutrality 
at equilibrium (the index '0' denotes the equilibrium 
quantities). Given this normalization, the ion viscos- 
ity rji is scaled into a (dimensionless) expression as 
V = Vi/i^pi^h)- As usual, m, denotes the ion mass, 
e is the electron charge, and ks is Boltzmann's con- 
stant; the ionic charge state was assumed to be Zi = 1. 
Ionic thermal effects are neglected, for simplicity. 

The electron distribution, assumed to be non- 
Maxwellian, is modelled by a re— type distribution 
[3TL [HI ; the electron density is thus given by: 



t<l> 



\)k B T e 



-re+1/2 



(4) 



where dimensions are reinstated at this step (only) for 
the (physical) electric potential $ and for the electron 
density h e (both latter quantities in physical dimen- 
sions, here, as opposed to the dimensionless counter- 
parts (j> and n). The (reduced) electron density in 
Poisson's equation ([3]) is thus taken to be of the form: 



-re+1/2 



(5) 



as results from the first momentum of the one- 
dimensional re-distribution function [34]]. Here, T e is a 
real constant, which corresponds to the electron tem- 
perature of the related Maxwellian plasma with the 
same particle and energy density (3~H . It is underlined 
that the re parameter holds a physical meaning if it 
belongs to the domain (3/2, 00); in the limit re — > 00, 
Eq. ([5]) retrieves the Maxwellian expression [34| . 
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The power law dependence depicted in Eq. ([5]) 
makes the differential equation of intractable an- 
alytically. In order to overcome this difficulty, we as- 
sume that any disturbance of the electric potential is 
small (i.e. <C ksTe/e). In this regime, a Taylor ex- 
pansion of Eq ([5]) around this parameter can be per- 
formed, allowing to express the electron density as 
n e ~ 1 + Q(f> + C2<p 2 + 0(</> 3 ). The expansion param- 
eters ci , C2 are only dependent upon k in the form: 



ci 



k- 1/2 
k- 3/2' 



C 2 = Cl 



1/2 



2{k- 3/2) 



(6) 



It is worth adding here that ci, c% >0Vk£ (3/2, 00). 
Truncating the expansion at the second order, Eq. ([3]) 
can be rewritten in the form: 



— -r w 1 - Hi + ci< 
oar 



(7) 



Since all the relevant quantities in the equations of 
interest now only refer to the ions, we will hereafter 
drop the subscript i for simplicity of notation. 



III. LINEAR DISPERSION RELATION 

A linear dispersion relation can be readily obtained 
by linearizing Eqs. $TJ, and ([7]), in the form: 



uj 2 = 



k 2 



k 2 + ci 

which reads, in physical units: 

k 2 ul 



k 2 + (^/Ac) 2 



(8) 



(9) 



Eq. ([9]) carries two fundamental consequences. The 
first is that u> is real for all k, implying that no 
damped mode is present. This is due to the fact that 
a weak value of the damping coefficient was assumed, 
by scaling the ion viscosity as 77 = e 1 ^ 2 ?7o (cf- the 
next Section) , hence the latter does not appear in the 



dispersion relation. Furthermore, the plasma screen- 
ing length of a K-distributed plasma is a function of 

k: \p = \/(k — 3/2)/ (k — 1/2) Ad, where A^> repre- 
sents the screening length of a Maxwcllian distributed 
plasma with same temperature and density (see Fig. 
QJ. The non-thermality of the plasma has thus the 
effect of reducing the plasma screening length: the 
Maxwellian value is in fact retrieved in the limit case 
ft — > CO. 

It is worth noticing here that the dispersion relation 
obtained here coincides with Eq. (12) in Ref. [35j . 
disregarding the magnetic field and setting v = c\ 
therein. 

In the long wavelength limit (i.e., for k <C 1), Eq. 
© reads 



k 



Cl c s 



J ph 



(10) 



where the k— dependent phase speed is defined as 
v ph = W c i c s ■ The sound speed therefore deviates 
from the ion sound speed in the related Maxwellian 
plasma (c s ) by a factor c x . This implies that the 
sound speed (and subsequently the phase speed of an 
ion-acoustic wave) is reduced in the presence of su- 
perthermal electrons. 



IV. DERIVATION OF THE K-DVB 
EQUATION 

We proceed now to derive a ^-dependent evolu- 
tion equation for propagating localized (constant pro- 
file) disturbances of the plasma dynamic variables, by 
adopting a multiscale perturbation technique. This 
technique operates by stretching the spatial and tem- 
poral coordinate with the help of an infinitesimal pa- 
rameter e 0, [ljj : 



,1/2 



(x - v ph t) 



f 3/2 



(11) 



where v p h is the ion excitation speed. The dependent 
variables n, u and cj> are perturbed around their equi- 
librium values in a power scries in e as: 
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FIG. 1: (Color online) Dependence of the normalized re- 
dependent plasma screening length with k. 



n=l + eni + e 2 n-2 + e 3 n^ 
u = em + e 2 u 2 + e 3 uz + ■ ■ 



• = e0i + e 2 



(12) 



We consider a weak damping by scaling the ion kine- 
matic viscosity as i] = e 1 / 2 ^. 

Substituting Eqs. (p} and CH]) into Eqs. ©, (0 
and (0, and separating the lowest order of e (i.e. 
terms proportional to e 3 / 2 ) the first order perturba- 
tions are obtained as: n\ = 4>i/v 2 h , ui = <pi/v p h and 

Vph = y/l/ci. The latter expression outlines the fact 
that, at first order, the ion-acoustic shock propagates 
at the K-dependent sound speed (as defined in Eq. 

m 
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(|T3|) will be mostly dictated by the interplay of the 
plasma dispersion and dissipation (depicted, in our 
model, by the parameters B, C). Moreover, given the 
exclusive relation of such parameters to k and to the 
ion viscosity 77, it is reasonable to anticipate a certain 
relation between these two plasma quantities, in order 
to distinguish different regimes. In the following Sec- 
tion, we will therefore integrate Eq. (fT3"|). and study 
the stability of the obtained solutions, in different sce- 
narios of arbitrary and strong dissipation. 



FIG. 2: (Color online) Dependence of the nonlinear coeffi- 
cient A and the dispersion coefficient B upon the param- 
eter K. 



Proceeding to the second order of perturbation (i.e. 
terms proportional to e 5 / 2 ) a KdVB equation is ob- 
tained in the form: 



9<t>i , . , 9<f>i 



D 



C- 



(13) 



where the nonlincarity coefficient A, dispersion coeffi- 
cient B and dissipation coefficient C are defined as: 



B 



2k-1 



2k-3 ' 
-3/2 



2k-3J 



(14) 



C — 2a 
2 ■ 



In the dissipationless case, i.e. in the limit t\q = 0, 
one can get the expected KdV equation and the cor- 
responding nonlinearity and dispersion coefficients as 
reported in Ref. for superthermal plasmas, upon 
substitution with a = Ti/T e = therein (in account 
of our cold ion fluid hypothesis). It is clear from 
the above expressions that both the nonlinear term 
A and the dispersion term B are strongly influenced 
by k. The nonlinear coefficient A is higher for stronger 
superthermality (lower «), while the dispersion co- 
efficient B is lower for stronger superthermality, as 
shown in Fig. [2] The expected Maxwellian limit 
A = 1 and B = 1/2 is recovered for k — > 00; this 
is identical to, e.g., Eq. (28) in Ref. [45| (consider- 
ing a = A, (3 = B, 5 = riio/n e o — 1 and a = 
therein). It is worth noting again that both A, B > 

Vk e (3/2,oo). 



V. SHOCK-PROFILE SOLUTIONS OF THE 
K-DVB EQUATION 

The nature of the electrostatic shock-like solutions 
of the Eq. (|13JI strictly depends upon the relation be- 
tween the equation coefficients A, B, C . In particular, 
given the overall weak dependence of A upon the pa- 
rameter k [see Eqs. (JTJ| an( l Fig. [2], it is reasonable 
to assume that the behavior of the solutions of Eq. 



A. General case: an analytical solution 

We consider here the general case in which neither 
dissipation nor dispersion can be neglected in solving 
the KdVB equation shown in Eq. (fT3|) . First of all, 
in order to ease the derivation, we will adopt a refer- 
ence frame moving with the shock; this is translated 
in transforming the spatial and temporal coordinates 
into: £ = ct(£ — Vt) and r = r (the time dependence 
disappears, since we anticipate stationary profile so- 
lutions). Here V represents the deviation of the shock 
speed from the plasma ion-acoustic speed and a" 1 
physically represents the shock width: qualitatively 
speaking, smaller values of a account for spatially ex- 
tended shock profiles, and vice versa. With this trans- 
formation of coordinates, Eq. (ff~3|) reads: 

Integrating once and imposing boundary conditions 
of the form: lim^oo </>i , dfyxjdC, , d 2 4>i/d( 2 = 0, Eq. 
([Pol) reduces to: 



Ba 



2 d 2 (/>i 
dC, 2 



d A x A 



0, (16) 



This equation is analytically solvable via the hyper- 
bolic tangent (tanh) approach [46Tj48j . The general 
solution takes the form: 



1 + tanh[(£ - Vt)/L ] 



(17) 



1 = 1QB/C represents the shock spa- 
D = $ max = 4$ n = I2C 2 /(25AB) 



where L = or 

tial width and u = v max = ito 
represents the shock amplitude (recall that AB > 
here). These values are obtainable by imposing that 
V = 6C 2 /(25B). This value of the velocity must 
be used in order to satisfy the boundary condition: 
lim^ +00 0i = V/A - 6C 2 /(25AB) = (Ifj It is 
worth noting here that this solution corresponds to a 
shock travelling towards positive £ (V > 0). It is read- 
ily seen that imposing opposite boundary conditions, 
a shock travelling in the opposite direction (V < 0) 
would be obtained with the same absolute values of 
width, amplitude and velocity. 
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A first interesting outcome of the analysis is that all 
the shock relevant parameters are functions of both re 
and r/ : 



1 



10 

no 

*0 = TM% 2 



2k- 
(2k- 



-3/2 



V 



100 ''0 (k-1)(2k-3) 
3/2 



(18) 



3_„2 / 2K-l \ 
% ^2^3 ) 



2o 



Recall that V is the shock velocity increment above 
the acoustic speed (i.e., the shock velocity is V s hock = 
c s + V), which is related to the viscosity ?]q. It is in- 
teresting to note that higher deviations from a pure 
Maxwellian behavior (i.e., smaller re) implies shocks 
that are narrower, faster and with a more ample po- 
tential jump (see Fig. [3] where a reference value of 
the viscosity 7; = 1 has been considered). Nonethe- 
less, for given re, the product $o^o ^ s constant, thus 
retaining a well-known feature of soliton solutions of 
the KdV equation. 

We now proceed to investigate the stability of a 
small perturbation (0i) around the exact solution 
given in Eq. (|17|) in the form: cf> = </>i(£, r)+e<j)\, where 
e is an infinitesimal parameter (e <C 1). Substituting 
</> in Eq. (|16p and linearizing with respect to <f>\, we 
obtain a differential equation for the perturbation that 
reads: 



Ba 



2 d2 4>i 
d( 2 



Ca- 



+ 4i 1 (A<i> 1 -V) = . (19) 



This equation admits a general solution proportional 
to e pl », where the parameter p obeys the characteristic 
polynomial: 



5a V ~ Cap + A<j) 1 -V = 

_ C±^C 2 -AB(A<t n -V) 



(20) 



J> = 



2Bc 



Perturbations around the exact solution will there- 
fore present an exponentially increasing (decreasing) 



behavior if p is positive (negative) or an oscillatory 
behavior if p is imaginary; in other words, the shock 
would sustain oscillatory perturbations if the discrim- 
inant of Eq. (j2T))) is negative: 



C 2 - 4B(A(j)i - V) < 



(21) 



Substituting the expression for 4>i given in Eq. (|17[) . 
such a condition is translated into: 



1 + 12[l + tanh(C)] 2 < 



(22) 



which is never satisfied, regardless of the value of 
The analysis therefore suggests that, in the general 
framework, oscillatory shock fronts are not allowed, 
regardless of the plasma characteristics. 

The two roots of Eq. (|20|) must be therefore real: 
P\-,P2 G K- The formal solution of Eq. (fH))) can thus 
be written as: 



F x e f 



F 2 e 



P2C 



(23) 
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FIG. 4: (Color online) Temporal evolution of the ion- 
acoustic shock described by Eq. (|17|l in a Maxwellian 
plasma (k — > 00). C = ^ = 0.5 and V = 0.1 have been 
considered. 
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FIG. 3: (Color online) Dependence of the shock width 
Lo (dashed red line), shock amplitude $ m „ = 4<E>o (dot- 
dashed blue line) and shock velocity V (solid green line) 
upon the non-thermality parameter k assuming r)o = 1. 



FIG. 5: (Color online) Temporal evolution of the ion- 
acoustic shock which is a solution for C = 2ja = 2 prop- 
agating in a superthermal plasma (« = 3), which is less 
dissipative (C = ^ = 0.2). 
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where F\ and F2 are integration constants. Recall- 
ing Eq. (f2Tj)) . it is readily obtainable that the sum of 
the two roots p\,P2 must be equal to Ca/(Ba 2 ) = 
IO770 > 0. The sum being positive, indicates that at 
least one of the roots must be positive implying that 
a small perturbation around the exact solution will be 
exponentially increasing in £, therefore disrupting the 
shock. No perturbations of the shock profile can thus 
be supported in this regime. In other words, the so- 
lution (|17[) provides a strictly monotonic profile (it is 
straightforward to show that the sign of the derivative 
of ipi is uniquely defined) which is unstable to external 
perturbations. 

A number of numerical simulations have been per- 
formed in order to verify our theoretical prediction. 
The propagation of the ion-acoustic shock has been 
analyzed for different plasma scenarios by numerical 
integration of the KdVB equation employing a Runge- 
Kutta 4 method. A time step of 10 -4 and a spatial 
grid with size 0.1 have been considered. As a test for 
the numerical code, the case of a shock propagating 
in a Maxwellian plasma (ft — > 00) and moderate dissi- 
pation 770 = 1 has been first studied numerically. The 
shock is in fact seen to propagate unperturbed (see 
Fig. [§. 

In order to trace the influence of the nonthermal- 
ity (kappa) and damping (770) parameter (s), we may 
consider a propagating shock crossing the interface be- 
tween two plasmas with either different k or different 
viscosity 770- in other words, we have considered the 
exact solution in plasma L, say, as initial condition 
to integrate the evolution equation (|T3"f in plasma R 
(here L and R, stand for the left and right plasma re- 
gions, respectively, assumed to differ in the value of 
one parameter only, either k or 770, all other parame- 
ters being the same). 

First, we have simulated the behavior of a shock 
crossing the interface between two plasmas with same 
k (here k = 3) but different viscosity (from rjo = 4 
to 770 = 0.4). The change in viscosity only induces 
a widening of the shock width, clearly seen in Figure 
which is brought about by the lower dissipation 
encountered by the shock. 

The situation is qualitatively different if we consider 
a transition in the superthermality of the plasma. If 
the shock penetrates into a region of lower k (i.e. from 
a Maxwellian to a k = 3 plasma, in our simulation) 
it is seen to preserve its strict monotonicity and its 
width; see Fig. Uta). On the other hand, the reverse 
phenomenon (i.e. a shock which is stable in a su- 
pcrthcrmal plasma entering in a Maxwellian environ- 
ment) presents a qualitatively different behavior. The 
shock in this case starts to develop a small hump in its 
spatial profile [see r = 250,400 in Fig. [f^b)]. These 
results can be interpreted in the following way. An 
increase in the n parameter (i.e. a smaller deviation 
from a pure Maxwellian behavior) implies an increase 
in the dispersive term B in the KdV-B equation [see 
Eqs. (TT4"1) ]. As it will be extensively discussed in the 



next section, oscillations at the front of an unstable 
shock occur only when the dispersion parameter dom- 
inates over the dissipative term. An abrupt increase 
of the dispersion parameter (as it is the case for an 
increase in the parameter k) can verify this condition, 
thus inducing oscillations at the shock front. On the 
other hand, the reverse case (transition from a Maxel- 
lian to a k = 3 plasma) carries the only consequence of 
further decreasing the dispersion parameter, thus fur- 
ther forbidding the onset of oscillations at the shock 
front. 

In the following, we proceed by adopting the same 
general treatment discussed in the previous section, 
yet considering the limit case of very weak (C 3> B) 
dispersion, compared to the damping coefficient C. 
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FIG. 6: (Color online) (a) Temporal evolution of an ion- 
acoustic shock (using as initial condition the stable so- 
lution (|17[) for a Maxwellian plasma), propagating in a 
strongly superthermal environment (ft = 3). (b) Temporal 
evolution of an ion-acoustic shock (using as initial condi- 
tion the solution which would be stable in a superthermal 
plasma, for n = 3), propagating in a Maxwellian environ- 
ment. Here C — ^ = 0.5 in both cases. 
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The opposite case (i.e. dominant dispersion) will not 
be treated here, since it will lead to a natural conver- 
gence to the well-known KdV equation. 

B. Weakly dispersive/strongly non-Maxwellian 
limit 

We now proceed by considering the particular case 
in which the dispersion of the background plasma 
(modelled by the parameter B) is negligible. Given 
the dependence of B upon k (Eq. (|14[) and Fig. [2]), 
this is translated into studying the behavior of shock- 
like solutions of the KdV-B equation in the case of 
strongly non-Maxwellian background plasmas [kappa 
in the vicinity of 3/2). In this particular regime, upon 
neglecting the dispersion coefficient B, Eq. (Til)]) be- 
comes: 

~ v ^c +A ^ = Ca W (24) 

which admits a general solution of the form (T^: 

MZ,t)=*iU- tanh[(C - Vt)/L^ . (25) 

Eq. ([25]) embodies a shock structure with speed V, 
amplitude $i = V/A and width L x = a' 1 = 2C/V. 
In order to check the stability of such solution a similar 
analysis as the one performed in the previous section 
will be performed. A perturbed solution of the form 
4> = 4>i (£, t) + e0i will be again considered and substi- 
tuted into Eq. <\'24\i after having performed a spatial 
integration with the same boundary conditions dis- 
cussed in the previous section. The linearization in 
<pi leads to a differential equation for the perturba- 
tion of the form: 

Ca^- - + V6 X = 0. (26) 

d( 

In this case, an analytical solution is obtainable by 
variable separation: 

0i=rsech 2 C, (27) 

where T represents the integration constant. Eq. 
(|2"T)l depicts a KdV soliton-type solution, suggesting 
that the limit case of dispersionless plasmas supports 
shocks with infinitcsimally perturbed fronts. If we 
return to an external fixed reference frame (charac- 
terized by the spatial coordinate £), the width and 
amplitude of the perturbation will be dictated by the 
interplay of C and V (remember that £ = (£ — Vt)/L\ 
with Li = 2C/V). In particular, higher propagation 
velocities (as much as smaller dissipation) will imply 
smaller perturbations. This conclusion is confirmed 
by numerical simulations of perturbed shock profiles 
shown in Fig. [7] for different values of ion viscosity 
[frame (a)] and propagation velocities [frame (b)]. 



It is interesting to test the behavior of shock solu- 
tions obtained in a purely dispersionless plasma, as 
they propagate through a region of not-negligible dis- 
persion. Again, given the dependence of B upon k 
(Eq. ([Ft]) and Fig. [2]), this refers to a scenario in 
which a shock propagating in a highly superthermal 
zone (k near 3/2) enters a "more Maxwellian" (higher 
k) region. Arguably, such a physical scenario reason- 
ably also models the temporal evolution of a shock 
through a plasma which is undergoing a progressive 
thcrmalization, which is a situation often encountered 
in experiments. Analytically, this can be tested by 
substituting Eq. ([25]) into Eq. (p~6]) . taking into ac- 
count that $i = V/A and L\ =2C/V. This leads to 
an equation of the form: 

BV 3 9 
- tanhC(tanh 2 C - 1) = . (28) 

This is equivalent to saying that a shock which is 
stable in the dispersion-less case, will preserve its sta- 
bility in a dispersive plasma only if the condition ([2"5]) 
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FIG. 7: (Color online) (a) Shock profile for different values 
of the dissipative term C, for a fixed velocity V = 1 in 
the regime of negligible dispersion, i.e. relying on solution 
(1251 1 with small perturbation; (b) Shock profile for different 
values of the velocity V, for fixed dissipation C = 0.1 in 
regime of negligible dispersion. Here, time r = 1 in both 
cases. 
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is satisfied. Although, strictly speaking, the latter 
condition is only satisfied (for arbitrary Q if B = 0, 
we may interpret the equality sign in the relation in 
a loose (~ 0) sense, thus allowing for different param- 
eter combinations which may approximately satisfy 
this condition, in practical terms. The most trivial 
cases to consider are either B = or C —> oo (which 
both recover the dispersionless evolution of the sys- 
tem) or tanh£ = ±1 (i.e. ( — > oo which suggests that 
a shock profile will keep to be stable in the regions far 
from the potential jump, regardless of the dispersion). 
However, a perturbation on the shock profile might be 
negligible also if BV 3 /2AC 2 <^ 1. Retrieving the de- 
pendence of A, B upon the non-thermality coefficient 
condition of the following form is obtained: 



C 



> 



k - 3/2 



V 3 / 2 a /2(k-1)(2k-1) 



(29) 



The above equation gives a K-dependent condition 
which determines the nature (monotonic or oscilla- 
tory) of a dispersionless shock front in a dispersive 
plasma (cf. Figs. l9lfT2"j) . 

In order to test the condition obtained above, a set 
of numerical simulations have been performed, and 
the results are shown in Figs. [S]-[T21 The simulations, 
which were obtained in the same setup as described 
above, aimed at following the evolution in time of a 
shock profile which is stable in a dispersionless plasma 
[i.e. the exact solution of Eq. ([24]) depicted in Eq. 
(|2"5j) ] when it is forced to propagate through a plasma 
with non-negligible dispersion [whose KdVB equation 
is shown in Eq. (|16[)]. As mentioned above, one would 
expect the shock to preserve its monotonicity as long 
as the dissipation parameter is dominant compared 
to the dispersion one. This statement can be easily 
tested by noticing that, in our model, the dispersion 
and the dissipation coefficients are solely functions 
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FIG. 8: (Color online) Threshold for the stability of 
dispersion-less shock fronts as a function of the non- 
thermality parameter k. The points depict the initial con- 
ditions for the numerical simulations shown in the follow- 
ing. 
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FIG. 9: (Color online) Temporal evolution of Eq. (fT6)l 
with initial condition given by Eq. (|25[l for a superthermal 
plasma (k = 3, upper panel) and a Maxwellian plasma 
(k —> oo, lower panel) for C = ^ — 0.6. The shock profile 
exhibits a purely monotonic structure in both cases. 



of k and of the ion viscosity parameter respectively 
[whose direct consequence is the threshold condition 
depicted in Eq. (|29[) ]. Different values of the viscosity 
r/o have been therefore tested (see the red points in 
Fig. [8]) in both approximately Maxwellian (k = 10) 
and strongly superthermal plasmas (n = 3) in order 
to lay above and below the critical condition for os- 
cillatory shocks. In all the simulations, a reference 
value of the velocity V — 1 has been assumed. For 
C = 0.6 (Fig. [5]) the shock is seen to be always mono- 
tonic, consistent with the analytical threshold which 
predicts the shocks being always monotonic regardless 
of k. Similarly, for C = 0.1 and C = 0.01 (Figs. [TO] 
and ITTj) the shock always presents an oscillatory be- 
havior, being below threshold for all values of k. In 
both these cases it is interesting to note that the case 
k = 3 comprises less pronounced oscillations than the 
Maxwellian case, being closer to the threshold curve 
depicted in Fig. [5] A peculiar situation can be found 
for C = 0.4. In this case the shock is expected to be 
monotonic for k = 3 (above threshold) and oscillatory 
for k — 10 (below threshold). Indeed this behavior is 
confirmed by the numerical results shown in Fig. 1121 
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FIG. 10: (Color online) Temporal evolution of Eq. (fl6j) 
with initial condition given by Eq. (125 p for a superthermal 
plasma (k = 3, upper panel) and a Maxwellian plasma 
(k — >• oo, lower panel) for C = ^ = 0.1. Oscillations are 
present in both cases. 



FIG. 11: (Color online) Temporal evolution of Eq. (Tl6|) 
with initial condition given by Eq. (|25f) for a superthermal 
plasma (k = 3, upper panel) and a Maxwellian plasma 
(k — s- oo, lower panel) for C = ^ = 0.01. In both cases, 
the shock presents large oscillations. 



VI. CONCLUSION 

We have presented a theoretical study of the prop- 
agation dynamics of ion-acoustic shock waves in a 
plasma (in a one-dimensional geometry) characterized 
by an excess in the superthermal electron population, 
where the electrons obey the k— distribution function, 
is presented. An unmagnetized electron-ion plasma 
has been assumed, with dissipation accounted for by 
considering the (constant) ion kinematic viscosity. 

A Korteweg-de Vries-Burgers equation was derived 
for the electrostatic potential via a multiscale pertur- 
bation technique, and then analyzed analytically and 
numerically. The nonlincarity and dispersion coeffi- 
cients are shown to be dependent upon the superther- 
mality index k, while the resulting damping coefficient 
is constant and only depends on the ion kinematic 
viscosity. Increasing the deviation from a Maxwellian 
behavior (i.e., considering smaller n values) simulta- 
neously induces a decrease of the dispersion coefficient 
and an increase of the nonlinearity counterpart. As a 
consequence, ion-acoustic shocks in non-thermal plas- 
mas are narrower, faster and of larger amplitude, as 
compared to the Maxwellian case. 

Different types of shock solutions were analytically 
shown to exist depending on the relation between the 



dispersion and the dissipation coefficients. Although 
a strictly monotonic profile solution is analytically ob- 
tained in the general cases, substantial differences be- 
tween the regimes of strong and weak dissipation have 
been found and confirmed by matching numerical in- 
tegrations of the Korteweg-de Vries-Burgers equation. 
The transition of ion-acoustic shocks through the in- 
terface between plasmas of different characteristics 
(different non-thermality or viscosity) has been nu- 
merically studied, showing a sensible deformation of 
the shock profiles. 

It appears appropriate here to add a comment to 
(formally) related work brought to our attention in the 
last stages of our study. Interestingly, a recent study 
has investigated superthermality effects on relativis- 
tic ion-acoustic shocks occurring in electron-positron- 
plasmas [50| . building upon a similar earlier investi- 
gation of Maxwellian plasmas [5l|. We feel like pay- 
ing due credit to Ref. admitting that our al- 
gebraic setting is partly recovered from the model 
adopted therein, in the appropriate (classical, vanish- 
ing positrons) limit. We stress nonetheless the physi- 
cal scope of our study at hand, which is distinct and 
presents no overlap with Ref. |5CI On the other hand, 
we find that a number of methodological glitches (and 
consequent erroneous results) orig inating in Ref. [5ll 
have been inherited by Ref. [50J, thus invalidating 
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FIG. 12: (Color online) Temporal evolution of Eq. (fT6)l 
with initial condition given by Eq. (|25|) for a superthermal 
plasma (k = 3, upper panel) and a Maxwellian plasma 
(« — s> cxd, lower panel) for C = % = 0.4. The shock 
presents oscillations only for the Maxwellian case. 



the (numerical, mainly) results therein. Therefore, 
not only the essential physics, but also the analytical 
formalism employed in those studies differ (s) signifi- 
cantly from our results and method herein. We have 
chosen not to burden the presentation here with an 
extensive reference to these latter references, as the 
details will be discussed elsewhere [48| . We note with 
interest that a later study by the same author [52j was 
dedicated to the effect of the (Tsallis) q— nonextensive 
distribution on the dynamics of ion-acoustic shocks. 
The Tsallis distribution that was employed in that 
study [4l| is not identical, nor amenable to the one 
resulting from the original Vasyliunas [29j kappa the- 
ory that we employ here (refer to the discussion in 
the Introduction). Furthermore, we regret to find that 
the forementioned analytical discrepancy [HI, [BtJ HJ 
arises in Ref. as well, invalidating the results of an 
otherwise interestingly conceived study. Concluding, 
neither the scope nor the details of our work present 
any overlap with the above studies. 

Our simplified theoretical model represents a small 
yet steady step towards the rigorous understanding 
of the behavior of ion-acoustic shocks in non-thermal 
environments, which appear to be of fundamental im- 
portance in a wide range of astrophysical and labo- 
ratory scenarios. A better adherence to realistic sce- 
narios of interest would require the relaxation of some 



assumptions adopted in the paper (such as taking into 
account of magnetic fields, multi-dimensionality and 
more realistic models for the ion viscosity) which will 
be subject of further work. 
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Appendix A: Derivation of the hybrid Korteweg 
-de Vries — Burgers (KdVB) equation 



Although the method employed here may strictly 
speaking not be original, a brief outline of the pertur- 
bativc methodology leading to the derivation of the 
KdVB equation ([I5|). 

Our starting point is the set of Eqs. ^ - An- 
ticipating a slow space, and a very slow time variation, 
according to the rationale underlying the KdV theory 
for acoustic electrostatic solitons @ , we introduce the 
coordinates £ = e x / 2 (x — v p ht) , 
writes 



e 3 / 2 i. One thus 



d_ 

dx 



,1/2 



d_ 



d_ 



+ e 



3/2_ 



dl 



(Al) 

Applying the stretching (|A1[) and the expansion 
(fT2|) into the model equations dU) - ([3j) , and isolating 
the contributions in e 3 / 2 , one obtains a set of equa- 
tions which can be solved to give the wave phase ve- 
locity 



Vph 



Cl 



where C\ was defined in the text. Interestingly, the lat- 
ter expression, here obtained as a compatibility con- 
dition, is exactly recovered via a different path if one 
considers the linear equations to obtain the ratio ui/k 
in the long wavelength limit; cf. the discussion in 
Section [IIII (This simply confirms the well-known su- 
peracoustic nature of electrostatic acoustic solitons.) 
The first-order density and velocity perturbation(s) 
are also obtained at this order, in terms of the electro- 
static potential (perturbation) 1; respectively given 
by 

m = <j>i/vl h , u\ = (f>i/v ph . 

In the next higher order in e (~ e 5 ^ 2 ), we obtain the 
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equations Now, differentiating Eq. (|A4[) with respect to £, one 

gets 

~ — 5 u p^-SF + + = ' ( A2 

9 ? u ph d Z d Z 

1 90i <9u 2 I d(f)i 902 % 9 2 0i 

+ + "W^^ ^ + ^-^^-20^ = 0. (A6) 



a2 JL 

-^|l + n 2 - Cl 2 - C 20? = . (A4) 

Multiplying Eq. (| A2|) by v p h then summing with Eq. 
l3|) in order to eliminate m 2 , one obtains 



Multiplying Eq. by ir^, adding to Eq. (fA3|) 

J_^l_ w 2 9»2 | 3 ^dcfrx | 90 2 7?q 9 2 0i = Q and then multiplying by Vph/2 leads us to write the re- 
Vph dr ph d£ Vp h <9£ d£ v p h <9£ 2 sultant equation in the form of the Korteweg-de Vries 

(A5) Burgers (KdVB) equation given in (|T3| . 
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